clear 

global state_max   statetaxcredit_t statemaxtaxcredit_t sunit_cap_t scc
%initialize parameters--------------------------------------------------
model_type='CollPercPolD'; % NoDem CollPercPolOwn 
%model_type='CollPercPolLL'
%model_type='CollPercPolTC'
net_metering=1;
%scc=41; % 41  148 330
spec="Opt"; % base specification for first set of tables




% model_type=strcat(model_type,string(scc));

%main talbes
all_specs=[ "OptCN"; "OptMD";  "OptTLCN" ;  "OptTLMD"; "Opt"; "Optmcf125"; "Optmcf150";    "OptNatCN"]; %   "OptTLCN" ; "OptTLMD"; "OptNS25";    "OptNatMD"; "Opt330";


%all_specs=[ "OptCN"; "OptMD";   "Opt"; "Optmcf125"; "Optmcf150"];

%line losses table and change model to LL
%model_type='CollPercPolDLL'
%all_specs=[ "OptCN"; "OptMD"; "Opt"];

%trans_con table and change model to TC
model_type='CollPercPolDTC'
all_specs=[ "OptCN"; "OptMD"; "Opt"];
%all_specs=[ "OptCN";  "Opt"];


%alt disc table and change model to CollPercPolDOD
%model_type='CollPercPolDOD'
%all_specs=[ "OptCN"; "OptMD"; "Opt"];


%new storage opt
%all_specs=[  "NROptNS25"; "OptNS25CN"; "OptNS25";  "NROptNS50";"OptNS50CN"; "OptNS50"; "NROptNS75";"OptNS75CN"; "OptNS75"];




%others
%all_specs=[ "OptCN"; "Opt"];

n_specs=size(all_specs,1);




npol=4; % types of pollutants, co2 nox pm25 so2
pol_types=["co2"; "nox"; "pm25"; "so2"];

[Sunit_state, Skwh_state, Scost_state, P_state, P_sale_state, PIns_state, PIns0_state, statemaxtaxcredit, statetaxcredit, sunit_cap, StateNames]=...
    read_state_data();



[A_t, M_t, BI_t, Nbar_i, HH_bin, CollPerc_t, Pol_t, Own_t, Inc_t, tract_nerc_xxwalk, tract_state_xxwalk,  distance_border, tract_border_xxwalk, population]=read_tract_data3();

ntract=size(A_t,1);
nbin=size(Nbar_i,2);


tract_state_xxwalk_mat=sparse(tract_state_xxwalk,1:size(tract_state_xxwalk,1),1)';

tract_state_xxwalk_mat=full(tract_state_xxwalk_mat); %ntract by nstate

statemaxtaxcredit_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*statemaxtaxcredit'),2)*ones(1,nbin);
statetaxcredit_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*statetaxcredit'),2)*ones(1,nbin);
sunit_cap_t=sum(tract_state_xxwalk_mat.*(ones(ntract,1)*sunit_cap'),2)*ones(1,nbin);



title='../Data/CleanedData/region_t.csv'; % census region xxwalk
tract_region_xxwalk=readmatrix(title);


load(sprintf('Outputs/Data%s%s', model_type, spec)) 

state_max=0;

%Make state level comparison table

Sunit_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_opt'),2);
Scost_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_opt'),2);
Skwh_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_opt'),2);

sub_amt_opt_t=calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
Sunit_opt_t, Scost_opt_t, Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t);  

sub_amt_opt_t=sub_amt_opt_t(:,1); 

Sunit_base_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_state'),2);
Scost_base_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_state'),2);
Skwh_base_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_state'),2);

state_max=1;

sub_amt_base_t=calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
Sunit_base_t, Scost_base_t, Skwh_base_t, P_t, PSale_t, PIns_t, PIns0_t);  

state_max=0;

sub_amt_base_t=sub_amt_base_t(:,1); % just one per tract and all bins are identical here

HH_state=sum((HH_t*ones(1,nstate)).*tract_state_xxwalk_mat)'; 



sub_base_state=sum(((HH_t.*sub_amt_base_t)*ones(1,nstate)).*tract_state_xxwalk_mat)'./(HH_state);
sub_opt_state=sum(((HH_t.*sub_amt_opt_t)*ones(1,nstate)).*tract_state_xxwalk_mat)'./(HH_state);


sub_base_nat=sum(HH_state.*sub_base_state)./sum(HH_state);
sub_opt_nat=sum(HH_t.*sub_amt_opt_t)./sum(HH_t);




BI_opt_state=sum((BI_opt*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
M_opt_state=sum((M_opt*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
Fisc_opt_state=sum(((Exp_sub_opt)*ones(1,nstate)).*tract_state_xxwalk_mat)'; 
D_opt_state=sum((D_opt*ones(1,nstate)).*tract_state_xxwalk_mat)'; 



diff_sub=100*(log(sub_opt_state)-log(sub_base_state));
diff_BI=100*(log(BI_opt_state)-log(BI_base_state));
diff_Fisc=100*(log(Fisc_opt_state)-log(Fisc_base_state));
diff_D=100*(log(D_opt_state)-log(D_base_state));


DiffTable=[StateNames diff_sub diff_BI diff_Fisc diff_D];
titles=["State" "diff_sub"  "diff_BI" "diff_Fisc" "diff_D"];

% append national ones
diff_sub_nat=100*(log(sub_opt_nat)-log(sub_base_nat));
diff_BI_nat=100*(log(sum(BI_opt_state))-log(sum(BI_base_state)));
diff_Fisc_nat=100*(log(sum(Fisc_opt_state))-log(sum(Fisc_base_state)));
diff_D_nat=100*(log(sum(D_opt_state))-log(sum(D_base_state)));

DiffTable=[DiffTable; "all" diff_sub_nat  diff_BI_nat diff_Fisc_nat diff_D_nat];

DiffTablesT=[titles; DiffTable];

title=(sprintf('Results/DiffTablesT%s%s.xls', model_type, spec));
writematrix(DiffTablesT, title); 



select_states=[4;  9; 19; 28; 33; 34; 45];

DiffTableTSelect=[titles; DiffTable(select_states,:)];

title=(sprintf('Results/DiffTableTSelect%s%s.xls', model_type, spec));
writematrix(DiffTableTSelect, title); 

subtable=[StateNames Skwh_opt Scost_opt Sunit_opt sub_opt_state]; 
titles=["State" "Skwh_opt" "Scost_opt" "Sunit_opt" "sub_opt_state"];

title=(sprintf('Results/subtable%s%s.xls', model_type, spec));
writematrix([titles; subtable], title);   
    

%levels talbe
LevelsTable=[StateNames sub_base_state sub_opt_state 1000*BI_base_state./HH_state 1000*BI_opt_state./HH_state];
titles=["State" "sub_base_state" "sub_opt_state" "BI_base_stateper1000HH" "BI_opt_stateper1000HH"];

%append national
LevelsTable=[LevelsTable; "all" sub_base_nat  sub_opt_nat 1000*sum(BI_base_state)./sum(HH_state) 1000*sum(BI_opt_state)./sum(HH_state)];

LevelsTableT=[titles; LevelsTable];
   
title=(sprintf('Results/LevelsTableT%s%s.xls', model_type, spec));
writematrix(LevelsTableT,title); 

LevelsTableTSelect=[titles; LevelsTable(select_states,:)];

title=(sprintf('Results/LevelsTableTSelect%s%s.xls', model_type, spec));
writematrix(LevelsTableTSelect,title ); 




%% amount of subsidy by type of subsidy

sub_opt_unit_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
Sunit_opt_t, 0*Scost_opt_t, 0*Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  

sub_opt_cost_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
0*Sunit_opt_t, Scost_opt_t, 0*Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  

sub_opt_kwh_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
0*Sunit_opt_t, 0*Scost_opt_t, Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  


state_max=1


sub_base_unit_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
Sunit_base_t, 0*Scost_base_t, 0*Skwh_base_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  

sub_base_cost_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
0*Sunit_base_t, Scost_base_t, 0*Skwh_base_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  

sub_base_kwh_nat=sum(HH_t.*calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
0*Sunit_base_t, 0*Scost_base_t, Skwh_base_t, P_t, PSale_t, PIns_t, PIns0_t))./sum(HH_t);  


state_max=0

Subtypes=[sub_base_unit_nat(1,1)./sub_base_nat sub_base_cost_nat(1,1)./sub_base_nat ...
    sub_base_kwh_nat(1,1)./sub_base_nat sub_base_nat./sub_base_nat; ...
   sub_opt_unit_nat(1,1)./sub_opt_nat sub_opt_cost_nat(1,1)./sub_opt_nat ...
   sub_opt_kwh_nat(1,1)./sub_opt_nat sub_opt_nat./sub_opt_nat ]'; %constant across all values of N_bar bins

Titles=["unit"; "cost"; "kwh"; "total"];

Subtypes=[Titles 100*Subtypes];

Titles=["" "Baseline" "Optimal"];


Subtypes=[Titles; Subtypes];

title=(sprintf('Results/SubsidyTypeNat%s%s.xls', model_type, spec));
writematrix(Subtypes,title ); 




%% saving stuff for big tables


DamagesBase=D_base;
DamagesPolBase=DPol_base;
%Tot_ut_sim=sum(Tot_ut_sim)
%WelfareBase=-D_base+Tot_ut_sim-sum(Exp_sub_sim);
UtBase=sum(Tot_ut_sim);
FiscalBase=sum(Exp_sub_sim);
BIBase=sum(BI_sim);
BIPHBase=sum(BI_sim)./sum(HH_t); % building per HH

DamagesOffBase=D0-D_base; % offset damages
DamagesOffPolBase=DPol0-DPol_base; % offset damages

WelfareBase=UtBase+DamagesOffBase-FiscalBase;

% by region
FiscalBaseR=zeros(4,1);
BIBaseR=zeros(4,1);  
BIPHBaseR=zeros(4,1);      
AveSubBaseR=zeros(4,1) ;   
NPanelsBaseR=zeros(4,1);
for r_i=1:4
   % FiscalOptR(r_i,1)=sum(Exp_sub_opt.*(tract_region_xxwalk==r_i));
    FiscalBaseR(r_i,1)=sum(Exp_sub_sim.*(tract_region_xxwalk==r_i));
    %BIOptR(r_i,1)=sum(BI_opt.*(tract_region_xxwalk==r_i));
    BIBaseR(r_i,1)=sum(BI_sim.*(tract_region_xxwalk==r_i));    
    
  %  BIPHOptR(r_i,1)=sum(BI_opt.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));
    BIPHBaseR(r_i,1)=sum(BI_sim.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));    
    AveSubBaseR(r_i,1)=sum(HH_t.*sub_amt_base_t.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));  

    NPanelsBaseR(r_i,1)=sum(M_sim.*(tract_region_xxwalk==r_i))./sum(BI_sim.*(tract_region_xxwalk==r_i));
end
    

%%  
%% alt specs
   
    UtOptSpec=zeros(n_specs,1); 
    DamagesOptSpec=zeros(n_specs,1);
    DamagesPolOptSpec=zeros(npol, n_specs);
    FiscalOptSpec=zeros(n_specs,1);
    BIOptSpec=zeros(n_specs,1);
    DamagesOffOptSpec=zeros(n_specs,1);
    DamagesOffPolOptSpec=zeros(npol, n_specs);
    BIPHOptSpec=zeros(n_specs,1);

    exceedSpec=zeros(n_specs,1);
    exceedSpec99=zeros(n_specs,1);
    
    FiscalOptSpecR=zeros(4, n_specs);
    BIOptSpecR    =zeros(4, n_specs);
    BIPHOptSpecR    =zeros(4, n_specs);
    AveSubSpecR =zeros(4, n_specs);
    NPanelsSpecR=zeros(4,n_specs); 
for c_i=1:n_specs
    
    spec=all_specs(c_i,1);
   
    load(sprintf('Outputs/Data%s%s', model_type, spec)) 
    
   % Tot_ut_opt=sum(Tot_ut_opt);

  %  WelfareOptSpec(c_i,1)=-D_opt+Tot_ut_opt-sum(Exp_sub_opt);
    UtOptSpec(c_i,1)=Tot_ut_opt;
    DamagesOptSpec(c_i,1)=D_opt;
    DamagesPolOptSpec(:, c_i)=DPol_opt;
    FiscalOptSpec(c_i,1)=sum(Exp_sub_opt);
    BIOptSpec(c_i,1)=sum(BI_opt);
    DamagesOffOptSpec(c_i,1)=D0-D_opt;
    DamagesOffPolOptSpec(:,c_i)=DPol0-DPol_opt;
    BIPHOptSpec(c_i,1)=sum(BI_opt)./sum(HH_t);

    exceedSpec(c_i,1)=meanExceed;
    exceedSpec99(c_i,1)=meanExceed99;

    Sunit_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_opt'),2);
    Scost_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_opt'),2);
    Skwh_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_opt'),2);

    if strcmp(spec,"OptTLCN")==1 || strcmp(spec,"OptTLMD")==1 
         load(sprintf('Outputs/TractLevel%s%s', model_type, spec)) 

        Sunit_opt_t= Sunit_tg;
        Scost_opt_t= Scost_tg;
        Skwh_opt_t= Skwh_tg;

    end

    sub_amt_opt_t=calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
    Sunit_opt_t, Scost_opt_t, Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t);  

    sub_amt_opt_t=sub_amt_opt_t(:,1); 



    for r_i=1:4
        FiscalOptSpecR(r_i,c_i)=sum(Exp_sub_opt.*(tract_region_xxwalk==r_i));
        BIOptSpecR(r_i,c_i)=sum(BI_opt.*(tract_region_xxwalk==r_i));   
        BIPHOptSpecR(r_i,c_i)=sum(BI_opt.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));
        AveSubSpecR(r_i,c_i)=sum(HH_t.*sub_amt_opt_t.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));  
        NPanelsSpecR(r_i,c_i)=sum(M_opt.*(tract_region_xxwalk==r_i))./sum(BI_opt.*(tract_region_xxwalk==r_i));
    end


end

WelfareOptSpec=UtOptSpec+DamagesOffOptSpec-FiscalOptSpec;

Counter=[BIBase FiscalBase DamagesBase DamagesOffBase; BIOptSpec FiscalOptSpec DamagesOptSpec DamagesOffOptSpec]; % ...
   %       BIOpt FiscalOpt DamagesOpt DamagesOffOpt];
%Counter=[BIBase FiscalBase DamagesBase DamagesOffBase; BIOptSpec FiscalOptSpec DamagesOptSpec DamagesOffOptSpec]; % ...   

titles=["Counterlist" "Installations" "Fiscal Cost" "Damages" "Damages Offset"];

names=["Baseline"; all_specs]; %; "Optimal Full"];
%names=all_specs';

Table=[names Counter];

Table=[titles; Table];



PV_nodep=20.5234564; % 25 years with 2% interest rate (no depreciation
% vertical table with regional stuf
CounterVert=round([1000*[BIPHBaseR BIPHOptSpecR; BIPHBase BIPHOptSpec' ]; zeros(1,1+n_specs); ...
    1*[AveSubBaseR AveSubSpecR ] ;zeros(1,1+n_specs); ...
  ... % 1/1000000*
    (1/PV_nodep)*1/1000*[FiscalBaseR FiscalOptSpecR ; FiscalBase FiscalOptSpec' ]; zeros(1,1+n_specs); ...
    (1/PV_nodep)*1/1000*[... %DamagesBase DamagesOptSpec' ;zeros(1,1+n_specs); ...
    DamagesOffPolBase DamagesOffPolOptSpec ; ...
    DamagesOffBase DamagesOffOptSpec' ; zeros(1,1+n_specs); ... 
    WelfareBase-WelfareBase WelfareOptSpec'-WelfareBase ; zeros(1,1+n_specs); ... 
    UtBase-UtBase UtOptSpec'-UtBase]  ;zeros(1,1+n_specs); ...
    NPanelsBaseR NPanelsSpecR],2);

rows=["Per1000HH Install:MW"; "Per1000HH  Install:NE"; "Per1000HH  Install:S"; "Per1000HH Install:W";  "Per1000HH Install:tot"; "zeros"; ...
    "Expected Sub:MW"; "Expected Su:NE"; "Expected Su:S"; "Expected Su:W";  "zeros"; ...
     "Tot Cost ($billions):MW"; "Tot  Cost($billions):NE"; "Tot  Cost($billions):S"; "Tot Cos($billions)t:W";  "Tot Cost($billions):tot"; "zeros"; ...
     ... ;"Damages($billions)"; "zeros";
     "annual co2 damages offset($millions)"; "annual nox damages offset($millions)"; "annual pm25 damages offset($millions)"; "annual so2 damages offset($millions)"; 
      "annual damages offset($millions)"; "zeros"; "annual welfare change ($millions)"; "zeros"; "annual HH utility ($millions)"; "zeros"; ...
     "Npanels Install:MW"; "Npanels  Install:NE"; "Npanels  Install:S"; "Npanels Install:W"];
 
 CounterVert=[rows CounterVert];
 
 names=[" " names'];
 
 CounterVert=[names; CounterVert];

 stop
 
 title=(sprintf('Results/CounterVert%s.xls', model_type));
%xlswrite(title, CounterVert); 
writematrix(CounterVert,title)
 


outof=ones(1, n_specs)*8750*7;
instances=exceedSpec'.*outof;

ExceedTable=[all_specs'; instances; outof; exceedSpec'; exceedSpec99' ]

 title=(sprintf('Results/ExceedTable%s.xls', model_type));
%xlswrite(title, CounterVert); 
writematrix(ExceedTable,title)


%%
%Alternative models
all_models=[ "NoDemD";   "CollPercD";  "CollPercPolD"; "CollPercPolOwnD"; "CollPercPolIncD"];
%all_models=[ "CollPercPolD";  "CollPercPolDIRH"; "CollPercPolDIRM"; "CollPercPolDIRL"; "CollPercPolDCC" ];

spec="OptCN";

%Note L M and H are opposite of paper (which refer to costs)
all_models=[  "CollPercPolD"; "CollPercPolDIRL"; "CollPercPolDIRM"; "CollPercPolDIRH"; "CollPercPolDCC"  ] %"CollPercPolCC"


n_model=size(all_models,1);

   % all of these are changes relative to baseline with current subsidies
   % but alternative utility func
    UtOptSpec=zeros(n_model,1); 
    DamagesOptSpec=zeros(n_model,1);
    DamagesPolOptSpec=zeros(npol, n_model);
    FiscalOptSpec=zeros(n_model,1);
    BIOptSpec=zeros(n_model,1);
    DamagesOffOptSpec=zeros(n_model,1);
    DamagesOffPolOptSpec=zeros(npol, n_model);
    BIPHOptSpec=zeros(n_model,1);
    AggBen=zeros(n_model,1);
    
    FiscalOptSpecR=zeros(4, n_model);
    BIOptSpecR    =zeros(4, n_model);
    BIPHOptSpecR    =zeros(4, n_model);
    AveSubSpecR =zeros(4, n_model);
    
for c_i=1:n_model
    
    model=all_models(c_i,1);
if scc>50
%    model=strcat(model,string(scc));
end
   
    load(sprintf('Outputs/Data%s%s', model, spec)) 
    
   % Tot_ut_opt=sum(Tot_ut_opt);

  %  WelfareOptSpec(c_i,1)=-D_opt+Tot_ut_opt-sum(Exp_sub_opt);
    UtOptSpec(c_i,1)=Tot_ut_opt-Tot_ut_sim;
    DamagesOptSpec(c_i,1)=D_opt-D_opt;
    DamagesPolOptSpec(:, c_i)=DPol_opt-DPol_sim;
    FiscalOptSpec(c_i,1)=sum(Exp_sub_opt)-sum(Exp_sub_sim);
    BIOptSpec(c_i,1)=sum(BI_opt)-sum(BI_sim);
    DamagesOffOptSpec(c_i,1)=D0-D_opt-(D0-D_sim);
    DamagesOffPolOptSpec(:,c_i)=DPol0-DPol_opt-(DPol0-DPol_sim);
    BIPHOptSpec(c_i,1)=sum(BI_opt)./sum(HH_t)-sum(BI_sim)./sum(HH_t);

    AggBen(c_i,1)=100*(D0-D_opt-(D0-D_sim))./(D0-D_sim);
    FiscalOptSpecPerc(c_i,1)=(sum(Exp_sub_opt)-sum(Exp_sub_sim))./sum(Exp_sub_sim);

    Sunit_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Sunit_opt'),2);
    Scost_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Scost_opt'),2);
    Skwh_opt_t= sum(tract_state_xxwalk_mat.*(ones(ntract,1)*Skwh_opt'),2);


    sub_amt_opt_t=calc_sub(N_mean*ones(size(Nbar_i)), A_t,  ...
    Sunit_opt_t, Scost_opt_t, Skwh_opt_t, P_t, PSale_t, PIns_t, PIns0_t);  

    sub_amt_opt_t=sub_amt_opt_t(:,1); 




    for r_i=1:4
     %   FiscalOptSpecR(r_i,c_i)=sum(Exp_sub_opt.*(tract_region_xxwalk==r_i));
     %   BIOptSpecR(r_i,c_i)=sum(BI_opt.*(tract_region_xxwalk==r_i));   
        BIPHOptSpecR(r_i,c_i)=sum(BI_opt.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i)) - ...
            sum(BI_sim.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));
        AveSubSpecR(r_i,c_i)=sum(HH_t.*sub_amt_opt_t.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i))- ...
        sum(HH_t.*sub_amt_base_t.*(tract_region_xxwalk==r_i))./sum(HH_t.*(tract_region_xxwalk==r_i));  
    end


end


%title=(sprintf('Results/nationalDamages%s.xls', model_type)); this table
%is not that useful
%writematrix(Table, title); 
%7	19.39792919
%20.52345647	19.39792919

%PV_nodep=20.5234564; % 25 years with 2% interest rate (no depreciation
% vertical table with regional stuf
CounterVert=round([1000*[ BIPHOptSpecR;  BIPHOptSpec' ]; zeros(1,n_model); ...
    AveSubSpecR  ;zeros(1,n_model); ...
      (1/PV_nodep)*1/1000*[... %DamagesBase DamagesOptSpec' ;zeros(1,1+n_specs); ...
     DamagesOffPolOptSpec ; ...
     DamagesOffOptSpec' ]; AggBen'; zeros(1,n_model); ...
     (1/PV_nodep)*1/1000*FiscalOptSpec'; 100*FiscalOptSpecPerc'],2);

rows=["Per1000HH Install:MW"; "Per1000HH  Install:NE"; "Per1000HH  Install:S"; "Per1000HH Install:W";  "Per1000HH Install:tot"; "zeros"; ...
    "Expected Sub:MW"; "Expected Su:NE"; "Expected Su:S"; "Expected Su:W";  "zeros"; ...
     "annual co2 damages offset($millions)"; "annual nox damages offset($millions)"; "annual pm25 damages offset($millions)"; "annual so2 damages offset($millions)"; 
      "annual damages offset($millions)"; "Change env ben percnet";"zeros"; "Annuitized Fiscal Cost ($millions)"; "% change gov cost"]; %; "zeros"; "annual welfare change ($millions)"; "zeros"; "annual HH utility ($millions)"];
 
 CounterVert=[rows CounterVert];
 
 names=[" " all_models'];
 
 CounterVert=[names; CounterVert];
 
 title=(sprintf('Results/CounterVertAltModels%s.xls', spec));
%xlswrite(title, CounterVert); 
writematrix(CounterVert,title)


